clc;
clear;

%% 迭代步骤 A 开始
% 读取网格数据
load msh

%% 迭代步骤 B 开始
% 设定边界条件
fy2 = 0; % 构建第一类边界条件数据 JB1
JB12 = [BP2, fy2*ones(size(BP2))]; % BP 及 BE 的定义见 2.3.2 小节
JB1 = [JB12];
dfy_dn1 = 0; % 构建第二类边界条件数据 JB2
dfy_dn3 = 0;
dfy_dn4 = 2;
JB21 = [BE1, ones(size(BE1(:,1)))*dfy_dn1, ones(size(BE1(:,1)))*dfy_dn1];
JB23 = [BE3, ones(size(BE3(:,1)))*dfy_dn3, ones(size(BE3(:,1)))*dfy_dn3];
JB24 = [BE4, ones(size(BE4(:,1)))*dfy_dn4, ones(size(BE4(:,1)))*dfy_dn4];
JB2 = [JB21; JB23; JB24];

clear JBV1 JBV2 JBV3 JBV4
clear BP1 BP2 BP3 BP4
clear JBP1 JBP2 JBP3 JBP4
clear P1 P2 P3 P4
clear BE1 BE2 BE3 BE4
clear JB12 JB21 JB22 JB23 JB24
clear thetax1 thetax2 thetax3 thetax4
clear thetay1 thetay2 thetay3 thetay4
clear dfy_dn1 dfy_dn2 dfy_dn3 dfy_dn4
clear fy2

%% 迭代步骤 C 开始
% 初始化总体系数矩阵 K 和 右边向量 F
K = zeros(N, N);
F = zeros(N, 1);

%% 迭代步骤 D 开始
% 逐个单元计算 Ke 并组合
for i = 1:E
    % 提取单元内三个节点的坐标, 存入 Jx 和 Jy
    for j = 1:3
        Jx(j) = JXY(JM(i,j),1);
        Jy(j) = JXY(JM(i,j),2);
    end
    % 计算三角形单元的面积
    Area = ((Jx(2)-Jx(1))*(Jy(3)-Jy(1)) - (Jx(3)-Jx(1))*(Jy(2)-Jy(1)))/2;
    % 计算单元中 a1,a2,a3,b1,b2,b3,c1,c2,c3
    a1 = (Jx(2)*Jy(3) - Jx(3)*Jy(2))/(2*Area);
    a2 = (Jx(3)*Jy(1) - Jx(1)*Jy(3))/(2*Area);
    a3 = (Jx(1)*Jy(2) - Jx(2)*Jy(1))/(2*Area);
    b1 = (Jy(2) - Jy(3))/(2*Area);
    b2 = (Jy(3) - Jy(1))/(2*Area);
    b3 = (Jy(1) - Jy(2))/(2*Area);
    c1 = (Jx(3) - Jx(2))/(2*Area);
    c2 = (Jx(1) - Jx(3))/(2*Area);
    c3 = (Jx(2) - Jx(1))/(2*Area);
    a = [a1, a2, a3];
    b = [b1, b2, b3];
    c = [c1, c2, c3];
    % 根据以上的 b,c 向量, 我们可以求出单元的单元系数矩阵
    for ki = 1:3
        for kj = 1:3
            Ke(ki, kj) = Area*(b(ki)*b(kj) + c(ki)*c(kj));
        end
    end
    % 总体方程组合
    for m = 1:3
        for n = 1:3
            K(JM(i, m), JM(i, n)) = K(JM(i, m), JM(i, n)) + Ke(m, n);
        end
    end
end

%% 迭代步骤 E 开始
% 逐个单元计算 Fe 并组合
for i = 1:length(JB2(:,1))
    II = JB2(i,1); % 提取单元序号
    boundary_side_number = JB2(i,2); % 提取边界边序号
    Fe = [0; 0; 0]; % Fe 初始化
    for j = 1:3 % 提取单元节点坐标
        Jx(j) = JXY(JM(II,j),1);
        Jy(j) = JXY(JM(II,j),2);
    end
    %% 单元边号为 1 时计算 Fe
    if boundary_side_number == 1
        q1 = JB2(i,5);
        q2 = JB2(i,5);
        L = sqrt((Jx(1)-Jx(2))^2 + (Jy(1)-Jy(2))^2);
        Fe(1,1) = L/6*(2*q1 + q2);
        Fe(2,1) = L/6*(q1 + 2*q2);
    end
    %% 单元边号为 2 时计算 Fe
    if boundary_side_number == 2
        q2 = JB2(i,5);
        q3 = JB2(i,5);
        L = sqrt((Jx(2)-Jx(3))^2 + (Jy(2)-Jy(3))^2);
        Fe(2,1) = L/6*(2*q2 + q3);
        Fe(3,1) = L/6*(q2 + 2*q3);
    end
    %% 单元边号为 3 时计算 Fe
    if boundary_side_number == 3
        q3 = JB2(i,5);
        q1 = JB2(i,5);
        L = sqrt((Jx(3)-Jx(1))^2 + (Jy(3)-Jy(1))^2);
        Fe(3,1) = L/6*(2*q3 + q1);
        Fe(1,1) = L/6*(q3 + 2*q1);
    end
    %% F 的组合
    for s = 1:3
        F(JM(II,s), 1) = F(JM(II,s), 1) + Fe(s,1);
    end
end

%% 迭代步骤 F 开始
% 采用乘大数代入法, 代入 JB1
dashu = 1e16;
for i = 1:length(JB1(:,1))
    K(JB1(i,1), JB1(i,1)) = K(JB1(i,1), JB1(i,1))*dashu;
    F(JB1(i,1), 1) = JB1(i,2)*K(JB1(i,1), JB1(i,1));
end

%% 迭代步骤 G 开始
% 求解方程, 显示结果
sudushi = K\F;
N = N;
E = E;
result = [JXY, sudushi]
JM

